Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 



Printed 7 December 2012 



(MN I^TeX style file v2.2) 



Solution to the Sigma-Problem of Pulsar Wind Nebulae. 



Oliver Porth 1 ' 2 *, Serguei S. Komissarov 1 ! , Rony Keppens 2 

1 Department of Applied Mathematics, The University of Leeds, Leeds, LS2 9GT 

2 Centre for mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Heverlee, Belgium 



Received / Accepted 



ABSTRACT 

We present first results of three dimensional relativistic magnetohydrodynamical sim- 
ulations of Pulsar Wind Nebulae. They show that the kink instability and magnetic 
dissipation inside these nebulae may be the key processes allowing to reconcile their 
observations with the theory of pulsar winds. In particular, the size of the termination 
shock, obtained in the simulations, agrees very well with the observations even for 
Poynting-dominated pulsar winds. Due to magnetic dissipation the total pressure in 
the simulated nebulae is particle-dominated and more or less uniform. While in the 
main body of the simulated nebulae the magnetic field becomes rather randomized, 
close to the termination shock, it is dominated by the regular toroidal field freshly- 
injected by the pulsar wind. This field is responsible for driving polar outflows and 
may explain the high polarization observed in pulsar wind nebulae. 
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1 INTRODUCTION 

It is now well established that Pulsar Wind Nebulae (PWN) 
are powered by relativistic winds from neutron stars, formed 
during violent deaths of massive normal stars. Theory of 
such winds, based on the physical processes in the magne- 
tospheres of neutron stars, predicts that they are strongly 
Poynting-dominated at their base (see Arons|2012[ and ref- 
erences therein). 

However, simple one-dimensional models of PWN fit 
the observations only if pulsar winds are particle-dominated. 
Quantitatively this conflict is best described using the wind 
magnetization parameter a defined as the ratio of the wind 
Poynting flux to its kinetic energy flux. Theoretical mod- 
els of pulsar magnetospheres and wind predict a ^ 1, 



that ideal MHD acceleration mechanisms can provide the 



whereas the ID models of PWN suggest a ~ 10 (Rees & 



Gunn||1974| |Kennel Coroniti|[T984l |Emmering fe Cheva- 
lier] 1987| |Begelman fe Li||1992| ). A slightly higher magneti- 
zation, a ~ 10 - ^, was later suggested by axisymmetric nu- 



merical simulations (Komissarov Sz Lyubarsky 2003 2004 



Del Zanna et al.|2004| |Bogovalov et al.|2005 ) 



Several possible explanations of this cr-problem have 
been put forward over the years. The simplest one is that 
the electromagnetic energy of the pulsar wind is converted 
into kinetic energy of the wind on its way from the pulsar to 
the wind termination shock, where freshly supplied plasma 
is injected into the PWN. Although claims have been made 
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required energy conversion (e.g. Vlahakis 



become clear that this is not the case (e.g. 
2009) |Lyubarsky|[2009l [2010] 



2004), it has now 
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The acceleration can be facilitated via non-ideal pro- 
cesses, involving magnetic dissipation, in pulsar winds 



( Coroniti 1990| . For the wind of the Crab pulsar, the dissipa- 
tion length scale still significantly exceeds the radius of the 
wind termination shock (|Lyubarsky &; Kirk 2001 ), unless the 



pulsar produces much more plasma compared to the predic- 
tions of the current models of pair-production in pulsar mag- 
netospheres (see Arons|2012 and references therein) . Alter- 
natively, the energy associated with the alternating compo- 
nent of magnetic field of the striped wind can be rapidly 
dissipated at the termination shock itself ([Lyubarsky 2003; 



Sironi &; Spitkovsky 2011| . Although the striped-wind model 



allows conversion of a large fraction of the total Poynting 
flux into the internal energy of PWN plasma, much more is 
needed to approach the target value of a ~ 10 -3 10 -2 . 
Indeed, the dissipation is confined to the striped zone of the 
wind and only the alternating component of magnetic field 
dissipates. Closer to the poles, the magnetization of pulsar 
wind plasma remains very high even after it crosses the ter- 
mination shock and, as a result, the overall magnetization of 
plasma injected into the nebula is much higher than that of 
the Kennel-Coroniti model, unless the striped zone spreads 
over almost the entire wind, implying that the pulsar is an al- 
most orthogonal rotator ( |Coroniti|p"9 90 ; Ko missarov|2012 ). 

In contrast to these ideas, Begelman (1998) argued that 
the ID models of PWN could be highly unrealistic and their 
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predictions not trustworthy. He has shown that the plasma 
configuration assumed in these models is unstable to the 
magnetic kink instability and speculated that the disrupted 
configuration may be less demanding on the magnetization 
of pulsar winds. Indeed, one would expect the magnetic pres- 
sure due to randomized magnetic field to dominate the mean 
Maxwell stress tensor, and the adiabatic compression to have 
the same effect on the magnetic pressure as on the thermo- 
dynamic pressure of relativistic gas. Under such conditions, 
the global dynamics of PWN produced by high-cr winds may 
not be that much different from those of PWN produced by 
particle-dominated winds. These expectations have received 



strong support from the recent numerical studies (Mizuno 



et al. 2011) of magnetic kink-instability of the cylindrical 



magnetostatic configuration, which was used in Begelman & 
[Li (1992) to model PWN. These simulations have shown a 
relaxation towards quasi-uniform total pressure distribution 
inside the computational domain on the dynamical time- 



scale. pVlizurioetaL (2011) also reported a significant mag- 



netic dissipation. 

The most important factor missing in the simulations by 



Mizuno et al. (2011) is the continuous injection of magnetic 



flux and energy in PWN by their pulsar winds. As a result, 
there is no termination shock whose size is an important 
parameter used to test theories of PWN against the obser- 
vational data. The observed strong polarization of PWN is 
suggestive of significant regular magnetic field in the form 
of concentric loops in their central regions, which could have 
been supplied by the pulsar wind relatively recently. The po- 
lar jets, most prominent in the Crab and Vela nebulae are 
best explained by the action of magnetic hoop stress of such 



fields ( Lyubarsky 2002 ; Komissarov &; Lyubarsky|2003 Del 



Zanna et al. 2004; Bogo valov et al.|2005| . Thus, the next log- 



ical step is to carry out 3D numerical simulations of PWN 
proper, similar in their setup to the previous axisymmetric 
simulations. We are currently conducting such studies and 
here report its first most interesting results. 



2 SIMULATION SETUP 

Initially, the computational domain is split in two zones sep- 
arated by a spherical boundary of radius r% = 10 18 cm. The 
outer zone is filled with radially-expanding cold supernova 
ejecta. The solution in this zone has Hubble's velocity profile 



(i) 



and constant density p e . The values of parameters p e 
and V[ are determined using the observed expansion rate 
of the Crab nebula and typical parameters of type-II super- 
novae. The inner zone is initially filled with unshocked pulsar 
wind. In our model of the wind, we assume that the alter- 
nating component of magnetic field in its striped zone has 
completely dissipated along the way from the pulsar to the 
PWN. Although this may not be the case and most of the 
dissipation occurs instead at the termination shock, dynam- 
ically this makes no difference (Lyubarsky 2003). The total 



energy flux density of the wind follows that of the monopole 
model flMichel|[l973t 




1000 



Figure 1. Electromagnetic luminosity over kinetic luminosity in- 
jected into the PWN against value before reconnect ion ctq for 
various angles of obliqueness a. If the obliqueness is larger than 
~ 45°, the luminosity fraction saturates at approximately this 
value. This is why our dynamical simulations can also provide 
reasonable models for highly magnetized cases with ctq > 10 3 . 
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The parameter b — 0.03 is added for numerical reasons. This 
energy is divided into magnetic, / m , and kinetic, /k, terms 



fm(r,6) = v(0) 



/tot(0,r) 



AM) 



/tot(r,0) 



(3) 



l + a{0)' JKy ' ' l + a(0) 

where a(9) is the latitude dependent wind magnetization. 

For simplicity, we assume that before dissipation of 
magnetic stripes the wind magnetization saturates at 

(6/8 ) 2 <t ; < O 
(To; > 0o 



d o {0) 



(4) 



where the additional (small) parameter Go — 10° was intro- 
duced to ensure that the Poynting flux vanishes on the axis 
for an axisymmetric model. Dissipation of magnetic stripes 
changes the magnetization of the striped wind zone so that 



(Tiff) 



*oi0)x a i0) 



l + *o(0)(l-x«(*)) 



(5) 



where 



(20«(0) I V • 



\ n/2-a<0 <n/2 + a 
otherwise 

(6) 

with <j> a (0) = arccos(— cot(0) cot (a)). Here, a signifies the 
magnetic inclination angle of the pulsar which determines 
the size of the striped wind zone (Komissarov 2012). Figure]!] 
shows the ratio of the wind total electromagnetic luminosity 
to its total kinetic luminosity, obtained via integration of the 
fluxes given by Eqs.|3] as a function of ao for different values 
of the magnetic inclination angle. In this letter we present 
simulations with ao = 0.01,1,3 and a = 45°. As one can 
see in Fig. [I] for ao = 3 the luminosity ratio is already very 
close to the asymptotic value in the limit ao — > oo, which 
is determined by the extent of the striped wind zone for all 
except very small values of a. This is why we expect our 
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Figure 2. Left panel: 3D rendering of the magnetic field structure in the model with cfq = 3 at t ~ 70 years after the start of the 
simulation. Magnetic field lines are integrated from sample points starting at r = 3 x 10 17 cm. Colors indicate the dominating field 
component, blue for toroidal and red for poloidal. Right panel: Azimuthally averag eda = (B^/B 2 )^ for the same model. 



models with gq — 3 not to be very different from those with 
much higher <7o, which can be found in real pulsar winds. 

The wind's magnetic field is purely azimuthal and 
changes direction at the equatorial plane. Its strength, as 
measured in the pulsar frame, is found via 

(7) 



B+ (r, 0) = ± vVMr, 0)/v , 



where v is the radial wind velocity. For simplicity, we adopt 
the same value of the wind Lorentz factor, r = 10, for all its 
streamlines. Hence the numerical wind density follows as 



P (r,e) = f k (r,e)/(r 2 c 2 v). 



(8) 



When the simulations begin, the discontinuity at n be- 
tween these two initial zones splits and the termination 
shock first rapidly moves inside towards the wind origin, re- 
flecting the artificial nature of our initial configuration. Soon 
after, it re-bounces and begins to expand at a much slower 
rate, gradually approaching its asymptotic self-similar po- 



sition ( |Rees &; Gunn 1974). For an unmagnetized wind, 
ao = 0, its equatorial radius at this stage evolves as 



ro 



1/2 



1 + 



V n t 



-1/2 



(9) 



where v n is the expansion rate of the PWN, and t is the time 
since the start of the simulations. One can see that the time- 
scale of transition to the self-similar expansion is ~ n/^i, 
which is ~ 210 yr for the parameters of our simulations. We 
do not yet fully reach the self-similar phase in 3D simulations 
and hence use this solution as a reference. 

The simulations are performed with MPI-AMRVAC 
(Keppens et al. 2012), solving the set of special relativis- 
tic magneto-hydrodynamic conservation laws in cartesian 
geometry. We employ a cubic domain with edge length of 
2 x 10 19 cm, large enough to contain today's Crab nebula. 
Outflow boundary conditions allow plasma to leave the sim- 
ulation box. The adaptive mesh refinement starts at a base 
level of 64 3 cells and is used to resolve the expanding neb- 
ula bubble with a cell-size of Ax = 1.95 x 10 16 cm (on the 



fifth level). Special action is taken to properly resolve the 
termination shock and the flow near the origin. To this end, 
additional grid levels centered on the termination shock are 
automatically activated, depending on the shock size. For 
the simulations shown here, the shock is thus resolved by 
3 — 4 extra grid-levels (resulting in 8 — 9 levels in total) on 
which we employ a more robust minmod reconstruction in 
combination with Lax-Friedrich flux splitting. 2D compar- 
ison simulations are performed with equivalent numerical 
setup and differ only in the use of cylindrical coordinates in 
the r, z-plane. 



3 RESULTS 



In basic agreement with the simulations by Mi zuno et al.| 
(2011), the highly organized coaxial configuration of mag- 
netic field, characteristic of previous 2D simulations of 
PWN, is largely destroyed in our 3D models. However, the 
azimuthal component is still dominant in the vicinity of 
the termination shock, in the region roughly corresponding 
to Crab's torus (see Figure [2]), which is filled mainly with 
"fresh" plasma which is just on its way from the termination 
shock to the main body of the nebula. As we have pointed 
out in the introduction, the emission of this plasma could 
be behind the strong polarization observed in the central 
region of the Crab nebula. This figure also shows predomi- 
nantly poloidal magnetic field in the outskirts of PWN, close 
to the equator. However, the magnetic field is rather weak 
in this region. 



Also in agreement with |Mizuno et ah (2011), the total 
pressure distribution of our 3D solutions is much more uni- 
form compared to 2D solutions of the same problem. As a 
result, the expansion of PWN in 3D is more or less isotropic, 
whereas in 2D the artificially enhanced axial compression 
due to the magnetic hoop stress promotes noticeably faster 
expansion in the polar direction. As one can see in Figure [3] 
by the end of the simulations, the 2D solution with ao = 3 
begins to exhibit a jet breakout, similar to those observed in 
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Figure 3. Left panel: The total pressure distribution in the 3D (left hemisphere) and the corresponding 2D (right hemisphere) solutions 
for the same model with ao = 3 at t ~ 70 years. Right panels: Velocity magnitude and direction for the same solutions, with the 2D one 
at the top and 3D one at the bottom. Note the scale difference between left/right panels. 



the earlier 2D simulations of highly-magnetized young PWN 



of magnetars (Bucciantini et al.||2007[ |2008|, with applica- 



tion to Gamma Ray Bursts. 

Figure [4] shows the evolution of the equatorial radius 
of the termination shock with time in all our simulations 
together with the analytical prediction from Eq. [9] We an- 
ticipated that in 3D the shock radius turns out similar to 
that given by the analytical model for unmagnetized wind, 
as this was suggested in Begelman ( 1998). This is indeed the 



case, but Figure [4] also shows that the size of the termina- 
tion shock exhibits only a weak dependence on the initial 
magnetization ao, once ao > 1. However, we also expected 
to find a much smaller shock radius in 2D simulations, as 
their symmetry prevents the development of the key process 
in the Begelman's theory, the kink instability. What we have 
actually found is that, although in our 2D simulations the 
shock radius is indeed smaller, the difference is not dramatic. 

An explanation for this result is suggested by Figure [5] 
which shows the ratio of magnetic to thermal energy of sim- 
ulated PWN in our numerical models. One can see that, 
not only do our 3D solutions exhibit significant magnetic 



dissipation, which agrees with Mizuno et al. (2011), but 



the 2D solutions do as well. Indeed, for 2D models with 
ao = 1,3 this parameter decreases from ~ 0.3, which is 



close to the mean value of freshly- injected plasma (Komis- 
2012| , to r>j 0.03. Such a low value means that the 



sarov 



total pressure in the nebula is dominated by the gas pres- 
sure, just like in models with particle-dominated winds, and 
this is why the shock radius is close to that of the unmag- 
netized model. The property of our 2D models which makes 
magnetic dissipation possible is the opposite orientation of 
magnetic loops in the northern and southern hemispheres, 
in contrast with most previous 2D simulations which were 
limited to only one hemisphere. Strong mixing between the 




0.05 



10 20 30 40 50 60 70 80 90 100 

t [years] 

Figure 4. Shock size in the equatorial plane for the two- (thin 
lines) and three-dimensional case (thick lines). The solid black 
curve labeled "H" follows the hydrodynamic prediction according 
to Eq. [9] 



two hemispheres, discovered already in | Camus et al. (2009), 



reduces the characteristic length scale of magnetic inhomo- 
geneities and speeds up the magnetic dissipation. The dis- 
placement of magnetic loops near the polar axis via the kink 
instability is an additional way of facilitating magnetic dis- 
sipation in 3D models. This agrees with the observed lower 
level of magnetic energy inside PWN in these models (see 
Figure [5| . 

The data presented in Figure[3]show that the axial com- 
pression remains a feature of our 3D-solution near the termi- 
nation shock and that the "tooth-paste" effect of hoop-stress 
due to freshly- injected azimuthal magnetic field of the pulsar 
wind is still capable of driving polar jets. But compared to 
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Figure 5. Ratio of the total magnetic energy to the total thermal 
energy of PWN as a function of time for 2D ( thin lines ) and 3D 
(thick lines) models. Horizontal dash-dotted lines indicate the ex- 
pected ratios at the self-similar phase in the absence of magnetic 
dissipation ( [Komissarov 2012). 



close to the termination shock and driving polar outflows, 
required to explain the Crab jet, and jets of other PWN. 
However, these are much more moderate than in 2D mod- 
els. It appears that the strong relativistic jets emerging in 
previous 2D simulations of magnet ar bubbles, with appli- 
cation to Gamma Ray Bursts, may well be artifacts of the 
imposed axial symmetry. 
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the 2D models, the region of operation of this jet-production 
mechanism is much more compact. 



4 DISCUSSION 

The results of our simulations provide strong support 
to Begelman's hypothesis that the essence of the sigma- 
problem of PWN relies in the inadequacy of using over- 
simplified models. He argued that these models were un- 
stable to kink modes and that randomization of magnetic 
field, resulting from this instability, would lead to a termi- 
nation shock of size similar to that of weakly-magnetized 
models, even if the actual magnetization was still very high. 
Our results confirm, and augment this by showing that the 
magnetization may actually be significantly reduced via on- 
set of magnetic dissipation accompanying randomization of 
magnetic field. The observed properties of PWN may mimic 
those of models with weakly magnetized winds not because 
magnetic stress of randomized magnetic field effectively re- 
duces to pressure but because the magnetic dissipation ren- 
ders the PWN gas-pressure-dominated (high-/3 plasma). In 
fact, high-magnetization of PWN plasma seems to be ruled 
out by the observations of synchrotron and inverse- Compt on 
emission of the Crab nebula, which show that the magnetic 
field is energetically sub-dominant to the population of rel- 
ativistic electrons by a factor of ~ 30 (Hillas 1998). 

Although the dissipation in the simulations ultimately 
is of numerical origin and occurs at the cell-size scale, the 
extreme resolution employed in a fully conservative, grid- 
adaptive computation can mimic processes creating such fine 
magnetic structures and hence the energy conversion has to 
be taken seriously. Moreover, the observed energy content of 
the Crab nebula does suggest efficient magnetic dissipation 
on a time scale which is significantly shorter than the age of 
the nebula but much longer compared to its Alfven crossing 
time (Komissarov 2012). 

The dynamics of the inner region of our 3D numerical 
models is still strongly influenced by the regular magnetic 
field of freshly-injected plasma. The hoop stress of this field 
is still capable of producing noticeable axial compression 
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